home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1996 January / macformat-033.iso / mac / Shareware City / Developers / VideoToolbox / VideoToolboxSources / Luminance.c < prev    next >
Encoding:
Text File  |  1995-09-19  |  39.5 KB  |  1,070 lines  |  [TEXT/CWIE]
  1. /*
  2. Luminance.c
  3. See Luminance.h for prototypes.
  4. See Luminance.note for explanation.
  5. Also see:
  6. D.G. Pelli and L. Zhang (1991) Accurate control of contrast on microcomputer displays. 
  7. Vision Research, 31:1337-1360.
  8.  
  9. Copyright (c) 1989-1995 Denis G. Pelli. This software is free; you may use it in
  10. your research and give it away to others, with the following restrictions. Any
  11. copy you give away must include this paragraph, unmodified, and, if you have
  12. made any changes to the rest of the file, must include a note, added to HISTORY
  13. below, giving your name, the date, and a description of the changes. This
  14. software may not be sold, whether in source or compiled form, without my
  15. permission. I hope you will find this software useful, but I can't promise that
  16. it will work for you, and am not offering any support. That's why it's free. I
  17. would appreciate reports of bugs and improvements.
  18.  
  19. Denis G. Pelli
  20. Department of Psychology
  21. New York University
  22. 6 Washington Place
  23. New York, NY 10003
  24. denis@cns.nyu.edu
  25.  
  26. HISTORY:
  27.  
  28. 4/24/89 dgp first working version.
  29. 4/28/89 dgp added device argument.
  30. 6/23/89 dgp added support for a video attenuator that has unequal gain in the three
  31.             channels. Added LoadClut() and LuminanceClut() routines.
  32. 7/30/89 dgp replaced call to SetEntries (which calls the Color Manager) by a call
  33.             to GDSetEntries (which calls the device driver).
  34. 8/10/89 dgp fixed bugs.
  35. 8/13/89 dgp rewrote LToV() to use nth order polynomial. Broke up LinearizeClut into
  36.             three routines--SetLuminanceRange, SetLuminance, and SetLuminances--the
  37.             last of which is equivalent to the old LinearizeClut.
  38. 8/16/89    dgp    After reading the Brooktree DAC manual I changed the documentation and
  39.             increased tolerance from 0.5 to 1 LSB, and now use the highest
  40.             luminance gain (worst case) over the requested range to transform it
  41.             into a luminance difference tolerance.
  42. 8/21/89    dgp    Made minor improvements in LToV() to deal with failure to converge due to 
  43.             getting stuck in a local minimum of the quartic. My solution is to now
  44.             require that Lmin and Lmax be filled in by the user in the calibration
  45.             structure. In practice this will guarantee that the minimum luminance
  46.             will not be below the local minimum found at the base of the rising curve.
  47.             Last week I also changed the LToV algorithm slightly, to do a bisection if
  48.             Newton's method would take us out of range.
  49. 9/7/89 dgp    Fixed bug in sorting routine, Sort().
  50. 9/10/89 dgp    Commented out the polynomial versions of VToL and LToV and replace them by
  51.             new versions that use a power law. The power law is a marginally better fit
  52.             than a 6th order polynomial (using only 4 instead of 7 parameters) and its
  53.             inverse can be computed much more quickly, about 0.3 ms instead of 2 ms.
  54. 9/11/89 dgp    Deciding to be indecisive, I introduced a conditional POWER_LAW_FIT
  55.             to control the choice of power law or polynomial fit for LToV and VToL.
  56. 9/16/89    dgp    Having finished writing the paper (Pelli & Zhang, 1991) I came back to
  57.             change the program to agree with what I wrote. I changed the equivalent
  58.             number tolerance to be the sum instead of the maximum of the 
  59.             variable-dac gains.
  60. 9/25/89 dgp Fixed minor bug in computation of tolerance. Now gives correct answer even
  61.             when the lowLuminance or highLuminance is out of range.
  62. 10/9/89 dgp    Made minor change so that when some of the gain[] are zero SetLuminance will
  63.             load zero into the unused colorSpec array elements. This is only a cosmetic
  64.             change.
  65. 10/28/89 dgp Made minor changes. I declared LToVPolynomial and LToVPower and
  66.             VToLPolynomial and VToLPower instead of having a conditional compilation. 
  67.             This makes both flavors of routine always available. I also changed 
  68.             LToVPolynomial to use LToVPower instead of LToVQuadratic for its initial 
  69.             guess. The quadratic guess was often poor and occasionally led to 
  70.             convergence problems. If the power law fit is not available then it reverts
  71.             to using the quadratic fit. If that's not available then it reverts to using
  72.             a middle-of-range guess.
  73. 10/28/89 dgp I eliminated the compile-time flag POWER_LAW_FIT and instead use whichever, 
  74.             polynomial or power, provides a better fit, as determined at run time. I 
  75.             biased the comparison of rms error to favor the power law fit since it has 
  76.             few parameters and can be inverted much more quickly. Note: if the 
  77.             power, polynomial, or quadratic fit parameters are not supplied then the 
  78.             appropriate LR.powerError, LR.polynomialError, or LR.quadraticError field 
  79.             should be set to infinity: INF or 1.0/0.0.
  80. 12/5/89    dgp    Added the trivial routine LToL() which enforces the bounds LP->LMin and 
  81.             LP->LMax.
  82. 4/2/90    dgp    Capitalized the word "To" in all function names, e.g. LToL().
  83. 4/22/90    dgp    Version 1.4. Made minor changes to the documentation. Renamed LToV() to
  84.             LToVFormulaic(), and introduced a new LToV() that uses a table to run faster.
  85.             LToV() takes an average of 170 µs, whereas LToVFormulaic() takes
  86.             1700 µs. I also now use the LuminanceTable, if available, to speed up
  87.             VToL() slightly, from 310 µs to 180 µs. Each call to SetLuminance() now
  88.             takes 0.9 ms whereas it used to take 2.5 ms. The loss in accuracy is
  89.             negligible. The interpolation error can be reduced still further by
  90.             increasing LUMINANCES_IN_TABLE. Each doubling of the table size quarters the
  91.             maximum possible error of the interpolation.
  92. 4/23/90    dgp    Added a few lines of code to LToV() to check near the last index, in the
  93.             spirit of Numerical Recipes in C, hunt.c, before starting the bisection
  94.             search. I updated the times in the 4/22/90 note to reflect the latest timing
  95.             by TestLuminance.
  96. 7/27/90    dgp    Tightened up the error checking for illegal entries into the clut. I have
  97.             the impression that, contrary to specification, the Apple video driver
  98.             crashes and corrupts itself if given an out-of-range entry value in a
  99.             setEntry Control call.
  100. 9/18/90    dgp    Changed all instances of "v" to "V". The final version of the Pelli
  101.             & Zhang (1991) paper refers to the nominal voltage v; this file
  102.             now refers to the "equivalent number" V; they are related by V=255*v.
  103. 9/24/90    dgp    Updated the documentation to correspond more closely to the (hopefully)
  104.             final version of the Pelli & Zhang (1991) paper.
  105. 10/29/90 dgp Doubled speed of SetLuminance(). Now SetLuminance() takes 133 ms to
  106.             do a complete clut for a small new luminance range, and takes 188 ms for a
  107.             large new range. This required minor changes to SetLuminanceRange() & LToV().
  108.             Changed all function headers to Standard C prototype style.
  109. 10/30/90 dgp Replaced phrase "clut value" by "nominal voltage", for consistency with
  110.             Pelli & Zhang (1991).
  111.             Introduced new fixed fraction data type called "Milli", defined in Milli.h.
  112.             By judicious replacement of double by Milli, particularly in the luminance
  113.             table, I have increased the speed of SetLuminance by a factor of three.
  114.             A Mac II now takes only 42 to 50 ms to build a whole clut. This new feature
  115.             is enabled by setting FAST_LUMINANCE to 1 in Luminance.h. There is a slight
  116.             increase in the luminance tolerance.
  117. 10/31/90 dgp More fine tuning. Now takes 31 to 44 ms to build a whole clut. Introduced
  118.             conditional FAST_LUMINANCE in Luminance.h that causes the same code to
  119.             be compiled with either Milli or double variables. The _SetLuminance()
  120.             routine is now quite polished, and runs fast. All variables and routines
  121.             whose names begin with underscore _ are written here with type Milli, but
  122.             if FAST_LUMINANCE is false then this type is redefined (in this file only)
  123.             as double, and all the Milli macros are redefined to be appropriate for
  124.             a double. So bear in mind that things beginning with underscore are of
  125.             unknown type. The virtue of being able to run the same code with Milli or
  126.             double computations is that if you suspect an error in the (homemade)
  127.             Milli arithmetic, you can try out double arithmetic simply by setting
  128.             FAST_LUMINANCE to 0 and recompiling.
  129.             In the interest of speed I have streamlined the tolerance
  130.             computation. The new tolerances are probably as good as the old ones.
  131. 11/2/90 dgp Added _SetLuminances() subroutine that substitutes linear interpolation
  132.             for most of the calls to _LToV(). This has greatly speeded up
  133.             SetLuminancesAndRange(), with minor loss of accuracy. A complete clut now
  134.             takes 9 to 31 ms. Even a Mac II can now compute low-contrast cluts
  135.             on the fly, between frames. (The user may wish to optimize the speed-
  136.             accuracy trade-off controlled by the LINEAR_V_DOMAIN parameter, which
  137.             determines the maximum width of the interpolation interval.)
  138. 11/6/90 dgp Replaced Milli by FIXED, which can be compiled as either double or Fixed.
  139.             The Fixed math is sufficiently precise as to give the same DAC values as
  140.             double math.
  141.             Now figure out optimal bit shift LT->L.LShift to preserve as much resolution
  142.             as possible when interpolating in _LToV(). 
  143.             A complete clut now takes 8 to 30 ms. For threshold stimuli the average
  144.             time will usually be less than 15 ms.
  145. 11/8/90 dgp    Eliminated gamma slope table since the speed-up it offered was too small
  146.             to measure. Tidied up the computation of LShift.
  147. 8/24/91    dgp    Made compatible with THINK C 5.0.
  148. 12/17/91 dgp Added new routines: IncrementLuminance() and GetV().
  149.             The former is used to obtain the lowest possible contrast. On a log contrast
  150.             scale, minus infinity (zero contrast) is a poor approximation to even
  151.             a small log contrast.
  152. 3/11/92    dgp    Minor editing of comments.
  153. 12/2/92 dgp Fixed minor bug in SetLuminanceRange() reported by Wei Xie and Ken Alexander
  154.             that could cause the THINK C Debugger to report an error when an 
  155.             uninitialized double V[] was converted to an int. The value was not used 
  156.             for anything, so this was innocuous.
  157. 12/17/92 dgp 1.94 Began enhancement to allow any dacSize, not just 8-bit dacs. I have
  158.             confirmed that everything still works fine with 8-bit dacs, but I don't
  159.             have a 9-bit dac video card to test it any further. •CAUTION: this probably
  160.             does not yet work correctly for 9-bit dacs. E.g. it seems likely that the 
  161.             fixed point math will break if VMax is increased to 511; a lot of careful 
  162.             error analysis went into the original coding, at which time I thought that 
  163.             it was safe to assume an 8-bit dac (i.e. VMax==255). •Removed obsolete 
  164.             support for THINK C 4. 
  165. 12/21/92 dgp No longer load unused dac bits.
  166. 1/4/93    dgp    LoadLuminances now checks the gdType.
  167. 5/10/93    dgp    If driver is in gray mode (i.e. not color) then LoadLuminances maps rgb to 
  168.             gray by simply copying the green component to the red and blue components.
  169. 5/12/93    dgp removed debugging code that forced isGray to always be true.
  170. 3/24/94    dgp    added LToE, LToEOrdered, and EToL.
  171. 6/29/94 dgp renamed "LtoEOrdered" to "LToEOrdered", fixing erroneous capitalization.
  172. 9/28/94    dgp    If necessary, IncrementLuminance now bumps up the upper luminance bound to 
  173.             include the incremented luminance.
  174. 10/25/94 dgp LoadLuminance and GetV now ignore a device with no driver. (E.g. a device
  175.             created by NewGWorld.)
  176. 12/29/94 dgp Fixed zero-divide halt in computing FAST_LUMINANCE _LToV returned value,
  177.             reported by Josh Solomon. Increased JMAX from 20 to 30 to eliminate spurious
  178.             warning in LToVPolynomial.
  179. */
  180. #include "VideoToolbox.h"        /* prototype for GDSetEntries() */
  181. #include "Luminance.h"
  182. static void Sort(int n,double arr[],int krr[]);    /* for internal use only */
  183.  
  184. #undef LongToFix
  185. #undef FixToLong
  186. #undef DoubleToFix
  187. #undef FixToDouble
  188. #if FAST_LUMINANCE
  189.     #define FIXED Fixed
  190.     #if !MACINTOSH
  191.         typedef long Fixed;
  192.         #define FixMul(x,y) DoubleToFix(FixToDouble(x)*FixToDouble(y))
  193.         #define FixDiv(x,y) DoubleToFix(FixToDouble(x)/FixToDouble(y))
  194.     #endif
  195.     #define LongToFix(x) ((long)(x)<<16)
  196.     #define FixToLong(x) ((x)>>16)
  197.     #define DoubleToFix(x) ((Fixed)((x)*65536.+0.5))
  198.     #define FixToDouble(x) ((double)(x)*(1./65536.))
  199.     #define D(x) FixToDouble(x)                            /* handy for debugging */
  200. #else
  201.     #define FIXED double
  202.     #define FixMul(x,y) ((double)(x)*(y))
  203.     #define FixDiv(x,y) ((double)(x)/(y))
  204.     #define LongToFix(x) ((double)(x))
  205.     #define FixToLong(x) ((long)(x))
  206.     #define DoubleToFix(x) ((double)(x))
  207.     #define FixToDouble(x) ((double)(x))
  208. #endif
  209.  
  210. #undef MAX
  211. #undef MIN
  212. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  213. #define MIN(a,b) ((a) < (b) ? (a) : (b))
  214.  
  215. double EToL(LuminanceRecord *LP,int entry)
  216. {
  217.     return GetLuminance(NULL,LP,entry);
  218. }
  219.  
  220. int LToE(LuminanceRecord *LP,double L,int firstEntry,int lastEntry)
  221. //Returns the index of the table entry in specified range with luminance closest to L.
  222. {
  223.     int i,entry;
  224.     double dL,latestDL;
  225.     
  226.     latestDL=INF;
  227.     entry=0;
  228.     for(i=firstEntry;i<=lastEntry;i++){
  229.         dL=fabs(L-GetLuminance(NULL,LP,i));
  230.         if(dL<latestDL){
  231.             entry=i;
  232.             latestDL=dL;
  233.         }
  234.     }
  235.     return entry;
  236. }
  237.  
  238. int LToEOrdered(LuminanceRecord *LP,double L,int firstEntry,int lastEntry)
  239. /*
  240. Returns the index of the table entry in specified range with luminance closest to L.
  241. Assumes ordered table from firstEntry to lastEntry, either increasing or decreasing.
  242. Returned entry is guaranteed to be in the range firstEntry...lastEntry.
  243. Based on Numerical Recipes "locate.c", page 117
  244. */
  245. {
  246.     register int lower,middle,upper;
  247.     register double LLower,LMiddle,LUpper;
  248.     int ascend;
  249.     
  250.     lower=firstEntry;
  251.     upper=lastEntry;
  252.     LLower=GetLuminance(NULL,LP,lower);
  253.     LUpper=GetLuminance(NULL,LP,upper);
  254.     ascend=(LUpper>LLower);
  255.     while(upper-lower>1){
  256.         middle=(upper+lower)/2;
  257.         LMiddle=GetLuminance(NULL,LP,middle);
  258.         if(L>LMiddle == ascend){
  259.             lower=middle;
  260.             LLower=LMiddle;
  261.         }
  262.         else{
  263.             upper=middle;
  264.             LUpper=LMiddle;
  265.         }
  266.     }
  267.     if(L-LLower<LUpper-L == ascend)return lower;
  268.     else return upper;
  269. }
  270. double SetLuminance(GDHandle device,LuminanceRecord *LP
  271.     ,int theEntry,double luminance
  272.     ,double lowLuminance,double highLuminance)
  273. /*
  274. Set one entry in the ColorSpec table (and the clut if device is not NULL) to
  275. a specified luminance. It's ok for lowLuminance to be greater than highLuminance; 
  276. they still designate a range.
  277. */
  278. {
  279.     FIXED _luminance;
  280.     
  281.     if(theEntry<0 || theEntry>=COLORS 
  282.         || device!=NULL && theEntry>=GDCLUTSIZE(device))
  283.             PrintfExit("\007SetLuminance: illegal entry %d\n",theEntry);
  284.     
  285.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  286.     _luminance=DoubleToFix(luminance);
  287.     _SetLuminance(LP,theEntry,_luminance);
  288.     if(device != NULL)LoadLuminances(device,LP,theEntry,theEntry);
  289.     return FixToDouble(_Tolerance(LP,_luminance));
  290. }
  291.  
  292. double SetLuminances(GDHandle device,LuminanceRecord *LP
  293.     ,int firstEntry,int lastEntry
  294.     ,double firstLuminance,double lastLuminance)
  295. /*
  296. Set a series of entries in the ColorSpec table (and the clut if device is not NULL)
  297. to a linear sequence of luminances. Assume this is the entire luminance range of 
  298. interest.
  299. */
  300. {
  301.     return SetLuminancesAndRange(device,LP,firstEntry,lastEntry
  302.         ,firstLuminance,lastLuminance,firstLuminance,lastLuminance);
  303. }
  304.  
  305. double SetLuminancesAndRange(GDHandle device,LuminanceRecord *LP
  306.     ,int firstEntry,int lastEntry
  307.     ,double firstLuminance,double lastLuminance
  308.     ,double lowLuminance,double highLuminance)
  309. /*
  310. Set a series of entries in the ColorSpec table (and the clut if device is not
  311. NULL) to a linear sequence of luminances. Uses last two arguments to set the
  312. luminance range.
  313.  
  314. I introduced a stack-space check before calling _SetLuminances(), since it calls
  315. itself recursively and may use 1000 bytes all together. However, printing an
  316. error message requires at least 4500 bytes. So I quit if the stack gets below
  317. 6000, so that the user will be presented with an understandable error message
  318. rather than a mysterious quit or crash.
  319. */
  320. {
  321.     FIXED _tolerance,_t,_dL64;
  322.     FIXED _firstL,_lastL,_firstV,_lastV;
  323.     
  324.     if(firstEntry<0 || firstEntry>lastEntry || lastEntry>=COLORS 
  325.         || device!=NULL && lastEntry>=GDCLUTSIZE(device) )
  326.         PrintfExit("\007SetLuminancesAndRange: illegal entries %d %d\n",firstEntry,lastEntry);
  327.     if(StackSpace()<6000)PrintfExit("SetLuminancesAndRange: not enough stack space.\007\n");
  328.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  329.     _firstL=DoubleToFix(firstLuminance);
  330.     _lastL=DoubleToFix(lastLuminance);
  331.     _firstV=_LToV(LP,_firstL + LP->_LOffset);
  332.     _lastV=_LToV(LP,_lastL + LP->_LOffset);
  333.     if(lastEntry != firstEntry){
  334.         #if FAST_LUMINANCE
  335.         _dL64=DoubleToFix(64.0*(lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  336.         #else
  337.         _dL64=DoubleToFix((lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  338.         #endif
  339.     }else _dL64=0;
  340.     _SetLuminances(LP,firstEntry,lastEntry,_firstL,_dL64,_firstV,_lastV);
  341.     if(device != NULL)LoadLuminances(device,LP,firstEntry,lastEntry);    
  342.     
  343.     /* quick and dirty estimate of tolerance */
  344.     _tolerance=LP->L._Lu[0] - _firstL;
  345.     _t=_lastL - LP->L._Lu[LUMINANCES_IN_TABLE-1];
  346.     if(_tolerance<_t)_tolerance=_t;
  347.     if(_tolerance<0)_tolerance=0;
  348.     _tolerance+=LP->_tolerance;
  349.     return FixToDouble(_tolerance);    /* estimate of max error in ΔL */
  350. }
  351.  
  352. void _SetLuminances(LuminanceRecord *LPtr,int first,int last
  353.     ,FIXED _firstL,FIXED _dL64,FIXED _firstV,FIXED _lastV)
  354. /*
  355. This routine eliminates most of the the slow _LToV() table lookups and
  356. interpolations by assuming that the gamma function is smooth enough that we may
  357. linearly interpolate over any interval of V no wider than LINEAR_V_DOMAIN. Since
  358. the gamma function is roughly parabolic, the consequent error in L will be
  359. proportional to the square of LINEAR_V_DOMAIN, and independent of where this
  360. interval is along V. Since we're inscribing straight line segments in a
  361. positively accelerating gamma function, we will overestimate luminance, and
  362. consequently V will be too low. The error will be greatest at the middle of each
  363. linearly interpolated V interval.
  364.  
  365. If the requested V interval is larger than LINEAR_V_DOMAIN, then it is split
  366. into two intervals by making a pair of recursive calls. LINEAR_V_DOMAIN is
  367. defined in Luminance.h. I suggest a value of 4, but it may be set larger to
  368. attain slightly higher speed at lower accuracy.
  369.  
  370. CAUTION: This routine is a stack hog. Each call uses up about 64 bytes on the
  371. stack, and it calls itself recursively, up to log2(256/LINEAR_V_DOMAIN) times.
  372. Do NOT change any of the declarations of the "new" variables used in the first
  373. if{} to "static", as that would cause the recursion to fail.
  374. */
  375. {
  376.     static int doLast=1;
  377.     int saveDoLast;
  378.     int newOne;
  379.     FIXED _newL,_newV;
  380.     register LuminanceRecord *LP=LPtr;
  381.     register int i,Vi;
  382.     register FIXED _VToGo,_g;
  383.     static RGBColor *RGBPtr;    /* static in order to minimize stack usage */
  384.     static int theEntry;        /* static in order to minimize stack usage */
  385.     static FIXED _V,_dV;        /* static in order to minimize stack usage */
  386.     register short leftShift;
  387.     
  388.     if(last-first>1 ){
  389.         _dV=_lastV-_firstV;
  390.         if(_dV>LongToFix(LINEAR_V_DOMAIN) || _dV<LongToFix(-LINEAR_V_DOMAIN)){
  391.             newOne=(first+last)>>1;
  392.             #if FAST_LUMINANCE
  393.                 _newL=_firstL+((newOne-first)*_dL64>>6);
  394.             #else
  395.                 _newL=_firstL+(newOne-first)*_dL64;
  396.             #endif
  397.             _newV=_LToV(LP,_newL + LP->_LOffset);
  398.             saveDoLast=doLast;
  399.             doLast=0;
  400.             _SetLuminances(LP,first,newOne,_firstL,_dL64,_firstV,_newV);
  401.             doLast=saveDoLast;
  402.             _SetLuminances(LP,newOne,last,_newL,_dL64,_newV,_lastV);
  403.             return;
  404.         }
  405.     }
  406.     if(last!=first)_dV=(_lastV-_firstV)/(last-first);
  407.     else _dV=0;
  408.     if(!doLast)last--;
  409.     leftShift=LP->leftShift;
  410.     for(theEntry=first,_V=_firstV;theEntry<=last;theEntry++,_V+=_dV){
  411.         /****** This section of code is copied from _SetLuminance() ****/
  412.         RGBPtr = &LP->table[theEntry].rgb;
  413.         *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  414.         _VToGo=_V - LP->_VFixed + LP->_VHalfStep;
  415.         for(i=LP->fixed;i<LP->dacs;i++) {
  416.             _g=LP->_gain[i];
  417.             Vi=_VToGo/_g;                        /* truncate to integer */
  418.             if(Vi>LP->VMax)Vi=LP->VMax;
  419.             if(Vi<LP->VMin)Vi=LP->VMin;
  420.             _VToGo -= Vi*_g;
  421.             ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  422.         }
  423.         /***************************************************************/
  424.     }
  425.     return;
  426. }
  427.  
  428. void _SetLuminance(LuminanceRecord *LPtr,int theEntry,FIXED _luminance)
  429. /*
  430. Set one entry in the ColorSpec table to a specified luminance. This is the
  431. private subroutine that actually does all the work for SetLuminance(). This
  432. routine is designed to run as fast as possible.
  433. */
  434. {
  435.     register LuminanceRecord *LP=LPtr;
  436.     register int i,Vi;
  437.     register FIXED _VToGo,_g;
  438.     RGBColor *RGBPtr;
  439.     register short leftShift=LP->leftShift;
  440.     
  441.     RGBPtr = &LP->table[theEntry].rgb;
  442.     *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  443.     _VToGo=_LToV(LP,_luminance + LP->_LOffset) - LP->_VFixed + LP->_VHalfStep;
  444.     for(i=LP->fixed;i<LP->dacs;i++) {
  445.         _g=LP->_gain[i];
  446.         Vi=_VToGo/_g;                        /* truncate to integer */
  447.         if(Vi>LP->VMax)Vi=LP->VMax;
  448.         if(Vi<LP->VMin)Vi=LP->VMin;
  449.         _VToGo -= Vi*_g;
  450.         ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  451.     }
  452. }
  453.  
  454. void IncrementLuminance(GDHandle device,register LuminanceRecord *LP,int theEntry)
  455. /*
  456. Make smallest possible increase of the luminance of one entry in the ColorSpec
  457. table. If necessary, bump up the upper luminance bound to include the incremented luminance.
  458. */
  459. {
  460.     double V,L,highL;
  461.     
  462.     V=GetV(device,LP,theEntry);
  463.     V+=LP->VHalfStep;
  464.     V+=LP->VHalfStep;
  465.     L=VToL(LP,V);
  466.     if(L>LP->highLuminance)highL=L;
  467.     else highL=LP->highLuminance;
  468.     SetLuminance(device,LP,theEntry,L,LP->lowLuminance,highL);
  469. }
  470.  
  471. FIXED _Tolerance(LuminanceRecord *LPtr,FIXED _luminance)
  472. {
  473.     register LuminanceRecord *LP=LPtr;
  474.     register FIXED _tolerance,_t;
  475.     
  476.     /* quick and dirty estimate of tolerance */
  477.     if(LP->L.exists==luminanceSet){
  478.         _tolerance=LP->L._Lu[0] - _luminance;
  479.         _t=_luminance - LP->L._Lu[LUMINANCES_IN_TABLE-1];
  480.         if(_tolerance<_t)_tolerance=_t;
  481.         if(_tolerance<0)_tolerance=0;
  482.         _tolerance+=LP->_tolerance;
  483.     }
  484.     else _tolerance=LP->_tolerance;
  485.     return _tolerance;
  486. }
  487.  
  488. double GetLuminance(GDHandle device,LuminanceRecord *LP,int theEntry)
  489. /*
  490. If device is not NULL then examines one entry in the actual clut, otherwise
  491. examines the ColorSpec table contained in *LP, and in either case returns the
  492. luminance that will be produced. Supplying an illegal entry value results in a
  493. returned value of -INF.
  494. */
  495. {
  496.     return VToL(LP,GetV(device,LP,theEntry));
  497. }
  498.  
  499. double GetV(GDHandle device,LuminanceRecord *LP,int theEntry)
  500. /*
  501. If device is not NULL (and has a driver) then examines one entry in the actual clut, otherwise
  502. examines the ColorSpec table contained in *LP, and in either case returns the V
  503. that will be produced. Supplying an illegal entry value results in a returned
  504. value of -INF.
  505. */
  506. {
  507.     RGBColor *myRGBPtr=NULL;
  508.     ColorSpec myCSpec;
  509.     double V;
  510.     int error;
  511.     
  512.     if(device != NULL && (*device)->gdRefNum!=0) {
  513.         if(theEntry<0 || theEntry>=GDCLUTSIZE(device)){
  514.             PrintfExit("GetLuminance: illegal entry %d\n",theEntry);
  515.         }
  516.         error=GDGetEntries(device,theEntry,0,&myCSpec);
  517.         if(error){
  518.             PrintfExit("GetLuminance: GDGetEntries error %d\007",error);
  519.         }
  520.         myRGBPtr=&myCSpec.rgb;
  521.     }else{
  522.         if(LP!=NULL && theEntry>=0 && theEntry<COLORS) 
  523.             myRGBPtr=&LP->table[theEntry].rgb;
  524.     }
  525.     if(myRGBPtr == NULL) return -INF;
  526.     V = LP->r*(myRGBPtr->red>>8);
  527.     V += LP->g*(myRGBPtr->green>>8);
  528.     V += LP->b*(myRGBPtr->blue>>8);
  529.     return V;
  530. }
  531.  
  532. void LoadLuminances(GDHandle device, LuminanceRecord *LP,
  533.     int firstEntry, int lastEntry)
  534. /*
  535. This just calls GDSetEntries() or GDDirectSetEntries() to load your ColorSpec
  536. table into the clut of your screen device. It is here simply to provide a
  537. cosmetic match to the call to SetLuminances(), for which loading the clut is
  538. optional. Note: if you prefer, instead of LP you may send just the address of a
  539. ColorSpec table, cast to (LuminanceRecord *), since a LuminanceRecord begins
  540. with a ColorSpec table.
  541. */
  542. {
  543.     int error;
  544.     short isGray,i;
  545.     RGBColor *RGBPtr;
  546.     ColorSpec table[256],*tablePtr;
  547.  
  548.     if(device==NULL || (*device)->gdRefNum==0)return;    // no device, or no device driver
  549.     isGray=!TestDeviceAttribute(device,gdDevType);
  550.     if(isGray){
  551.         for(i=firstEntry;i<=lastEntry;i++){
  552.             RGBPtr=&table[i].rgb;
  553.             RGBPtr->red=RGBPtr->green=RGBPtr->blue=LP->table[i].rgb.green;
  554.         }
  555.         tablePtr=&table[firstEntry];
  556.     }else tablePtr=&LP->table[firstEntry];
  557.  
  558.     switch((*device)->gdType){
  559.     case fixedType:
  560.         printf("LoadLuminances: this device has a fixed CLUT.\n");
  561.         return;
  562.     case clutType:
  563.         error=GDSetEntries(device,firstEntry,lastEntry-firstEntry,tablePtr);
  564.         break;
  565.     case directType:
  566.         error=GDDirectSetEntries(device,firstEntry,lastEntry-firstEntry,tablePtr);
  567.         break;
  568.     }
  569.     if(error) printf("\007LoadLuminances: error %d\n",error);
  570. }
  571.  
  572. double LToL(LuminanceRecord *LP,double L)
  573. {
  574.     if(L > LP->LMax) return LP->LMax;
  575.     if(L < LP->LMin) return LP->LMin;
  576.     return L;
  577. }
  578.  
  579. double VToL(LuminanceRecord *LP,double V)
  580. /*
  581. Return the luminance that would result from a nominal voltage.
  582. */
  583. {
  584.     return FixToDouble(_VToL(LP,DoubleToFix(V)));
  585. }
  586.  
  587. FIXED _VToL(LuminanceRecord *LP,FIXED _V)
  588. /*
  589. Return the luminance that would result from a nominal voltage.
  590. */
  591. {
  592.     register int i;
  593.     register LuminanceTable *LT;
  594.     register FIXED _di;
  595.     double LF,VF;
  596.     
  597.     LT=&LP->L;
  598.     if(LT->exists==luminanceSet){
  599.         _di=FixDiv(_V-LT->_VMin,LT->_dV);
  600.         i=FixToLong(_di);
  601.         _di -= LongToFix(i);
  602.         if(i<0)return LT->_Lu[0];
  603.         if(i>=LUMINANCES_IN_TABLE-1)return LT->_Lu[LUMINANCES_IN_TABLE-1];
  604.         return LT->_Lu[i]+FixMul(_di,LT->_Lu[i+1]-LT->_Lu[i]);
  605.     }
  606.     VF=FixToDouble(_V);
  607.     if(LP->powerError < 2.0*LP->polynomialError) LF=VToLPower(LP,VF);
  608.     else LF=VToLPolynomial(LP,VF);
  609.     return DoubleToFix(LF);
  610. }
  611.  
  612. double LToV(LuminanceRecord *LP,double L)
  613. {
  614.     return FixToDouble(_LToV(LP,DoubleToFix(L)));
  615. }
  616.  
  617. FIXED _LToV(LuminanceRecord *LP,FIXED _Lu)
  618. /*
  619. New, faster version uses a tabulated version of the gamma function VToL(). This
  620. replaces the old function LToVFormulaic(), formerly called LToV(). A subtlety is
  621. that I tabulate the gamma function _Lu[]=VToL() rather than its inverse, because
  622. this yields a much lower upper bound on the error of the linearly interpolated
  623. L. This is true because the gamma function is roughly parabolic. Equally spaced
  624. samples of a parabolic function yield equal maximum error of linear
  625. interpolation in all intervals. The cost of tabulating _Lu[] rather than its
  626. inverse is that we must search the table in order to find the largest _Lu[] less
  627. than or equal to L. The search time is minimized by using a search procedure
  628. that begins the search at the place found by the last search, since successive
  629. calls to LToV() are likely to request similar luminances.
  630. */
  631. {
  632.     register LuminanceTable *LT;
  633.     register int lo,hi,i;
  634.     register FIXED _dL,_dLTemp;
  635.     FIXED _V,_dLLo;
  636.     
  637.     LT=&LP->L;
  638.     if(LT->exists != luminanceSet){
  639.         /* get domain of the monotonic section of the gamma function */
  640.         LT->_VMin=DoubleToFix(LToVFormulaic(LP,LP->LMin));
  641.         LT->_VMax=DoubleToFix(LToVFormulaic(LP,LP->LMax));
  642.         LT->_dV=(LT->_VMax-LT->_VMin)/(LUMINANCES_IN_TABLE-1);
  643.         _V=LT->_VMin;
  644.         for(i=0;i<LUMINANCES_IN_TABLE;i++,_V+=LT->_dV) LT->_Lu[i]=_VToL(LP,_V);
  645.         /* just in case, impose monotonicity, from center on out */
  646.         for(i=LUMINANCES_IN_TABLE/2;i>0;i--)
  647.             if(LT->_Lu[i-1] > LT->_Lu[i]) LT->_Lu[i-1]=LT->_Lu[i];
  648.         for(i=1+LUMINANCES_IN_TABLE/2;i<LUMINANCES_IN_TABLE;i++)
  649.             if(LT->_Lu[i-1] > LT->_Lu[i]) LT->_Lu[i]=LT->_Lu[i-1];
  650.         #if FAST_LUMINANCE    /* compute the slope dL/dV */
  651.             /* find shift that maximizes ΔL precision without overflow */
  652.             _dL=0;
  653.             for(i=0;i<LUMINANCES_IN_TABLE-1;i++){
  654.                 _dLTemp=LT->_Lu[i+1]-LT->_Lu[i];
  655.                 if(_dL<_dLTemp)_dL=_dLTemp;
  656.             }
  657.             LT->LShift=Log2L(LONG_MAX/_dL);
  658.         #endif                
  659.         LT->exists=luminanceSet;
  660.         LT->latestIndex=-2;    /* invalid, so hunt will start from scratch */
  661.     }
  662.     /* hunt for L in LT->_Lu[] */
  663.     /* first check at latestIndex, latestIndex±1, ... */
  664.     lo=-1;
  665.     hi=LUMINANCES_IN_TABLE;
  666.     i=LT->latestIndex;
  667.     if(i>lo && i<hi){
  668.         if(_Lu>LT->_Lu[i])
  669.             {lo=i;i++;}
  670.         else
  671.             {hi=i;i--;}
  672.     }
  673.     /* simple bisection search, see Numerical Recipes in C, pages 98-99. */
  674.     while(hi-lo>1){
  675.         i=(hi+lo)>>1;
  676.         if(_Lu>LT->_Lu[i])lo=i;
  677.         else hi=i;
  678.     }
  679.     LT->latestIndex=lo;
  680.     if(lo<0){
  681.         LT->latestIndex=0;            /* nearest legal index */
  682.         return LT->_VMin;
  683.     }
  684.     if(lo>=LUMINANCES_IN_TABLE-1)return LT->_VMax;
  685.     _V=LT->_VMin+lo*LT->_dV;
  686.     _dL=_Lu-LT->_Lu[lo];
  687.     _dLLo=LT->_Lu[lo+1]-LT->_Lu[lo];
  688.     if(FAST_LUMINANCE){
  689.         if(_dL<LONG_MAX/LT->_dV)
  690.             _V+=LT->_dV*_dL/_dLLo;
  691.         else if((_dLLo<<LT->LShift) >= LT->_dV)
  692.             _V+=(_dL<<LT->LShift)/((_dLLo<<LT->LShift)/LT->_dV);
  693.     }else _V+=LT->_dV*(_dL/_dLLo);
  694.     return _V;
  695. }
  696.  
  697. double LToVFormulaic(LuminanceRecord *LP,double L)
  698. /*
  699. Find a nominal voltage that would give luminance L. Failing that, it returns the
  700. closest possible value. The answer is computed by inverting the formula for the
  701. gamma function.
  702. */
  703. {
  704.     if(LP->powerError < 2.0*LP->polynomialError) return LToVPower(LP,L);
  705.     else return LToVPolynomial(LP,L);
  706. }
  707.  
  708. double VToLPolynomial(LuminanceRecord *LP,double V)
  709. /*
  710. Return the luminance that would result from a nominal voltage. Uses a polynomial
  711. equation L = p[0] + p[1]*V^1 +... of any degree. dgp 8/11/89
  712. */
  713. {
  714.     register double L,VV;
  715.     register int i,m;
  716.     
  717.     m=LP->coefficients;
  718.     L=0.0;
  719.     VV=1.0;
  720.     for(i=0;i<m;i++){
  721.         L+=LP->p[i]*VV;
  722.         VV*=V;
  723.     }
  724.     return L;
  725. }
  726.  
  727. double VToLPower(LuminanceRecord *LP,double V)
  728. /*
  729. Return the luminance that would result from a nominal voltage.
  730. Uses a rectified power law:
  731.  
  732.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  733.     
  734.     Rectify(x)=x if x≥0
  735.     Rectify(x)=0 if x<0
  736.  
  737. dgp 10/28/89
  738. */
  739. {
  740.     double L;
  741.     
  742.     L=LP->power[1]+LP->power[2]*V;
  743.     if(L>0.0) L=LP->power[0]+pow(L,LP->power[3]);
  744.     else L=LP->power[0];
  745.     return L;
  746. }
  747.  
  748. double LToVPolynomial(LuminanceRecord *LP,double L)
  749. /*
  750. Find a nominal voltage that would give luminance L. Failing that, it returns the
  751. closest possible value. Solves a polynomial equation L = p[0] + p[1]*V^1 +... of
  752. any degree. First we find a quick approximate answer by solving the power law
  753. fit, LToVPower(). Then we use Newton's method to home in quickly on the solution
  754. to the polynomial fit. Newton's method is taken from Numerical Recipes in C. For
  755. a polynomial of degree 8 (which is what I recommend) this routine now does about
  756. 3? iterations and takes about 1.1? ms. The error in V is less than TOLERANCE,
  757. i.e. relative to full scale of 256 we get an accuracy of 14.6 bits. You can
  758. increase the TOLERANCE to reduce the number of iterations to make it slightly
  759. faster. A TOLERANCE of 1.0 (for 8-bit accuracy) takes 0.8? ms, not much of a
  760. savings in time, and a large loss in accuracy. Reducing TOLERANCE to 1e-8 yields
  761. an accuracy of nearly 35 bits and takes 1.5? ms, only slightly longer.
  762.  
  763. This routine will be called essentially once per clut entry, so computing a
  764. whole new clut with entries 0..255 will take 256*1.1? ms = 282? ms. Some
  765. experiments use only part of the clut, and will therefore take less time. In
  766. some of the experiments that do use the whole clut it may be desirable to
  767. compute the clut table in advance.
  768. dgp 8/11/89
  769. */
  770. {
  771.     #define TOLERANCE 0.01                    /* desired accuracy of V, see above */
  772.     #define JMAX 30
  773.     double a[MAX_COEFFICIENTS];
  774.     register int i;
  775.     int m,j;
  776.     register double f,df,u,VV,V,dV;
  777.     
  778.     if(L>LP->LMax) return LP->VMax;
  779.     if(L<LP->LMin) return LP->VMin;
  780.  
  781.     m=LP->coefficients;
  782.     if(m>MAX_COEFFICIENTS || m<1)
  783.         PrintfExit("LToVPolynomial: %d coefficients is too many or too few\007\n",m);
  784.     for(i=0;i<m;i++) a[i]=LP->p[i];
  785.     a[0]-=L;
  786.     if(LP->powerError < 0.2*LP->LMax) V=LToVPower(LP,L);    /* very good initial guess */
  787.     else {
  788.         if(LP->quadraticError < 0.2*LP->LMax) V=LToVQuadratic(LP,L);/* fair initial guess */
  789.         else PrintfExit("LToVPolynomial: neither power nor quadratic fit is available\007\n");
  790.     }
  791.     for(j=0;j<JMAX;j++){
  792.         /* evaluate function, i.e. the polynomial, and its derivative */
  793.         f=a[0];
  794.         df=0.0;
  795.         VV=1.0;
  796.         for(i=1;i<m;i++){
  797.             u=a[i]*VV;
  798.             df+=i*u;
  799.             f+=V*u;
  800.             VV*=V;
  801.         }
  802.         dV=f/df;    /* apply Newton's method */
  803.         if(V-dV<LP->VMin) dV=(V-LP->VMin)/2.0;    /* bisect */
  804.         if(V-dV>LP->VMax) dV=(V-LP->VMax)/2.0;    /* bisect */
  805.         V -= dV;
  806.         if(fabs(dV) < TOLERANCE) return V;
  807.     }
  808.     printf("LToVPolynomial(%7.2f): Warning, too many iterations. VToL(%5.1f)=%7.2f\n",V,L,VToL(LP,V));
  809.     return V;
  810.     #undef TOLERANCE
  811.     #undef JMAX
  812. }
  813.  
  814. double LToVPower(LuminanceRecord *LP,double L)
  815. /*
  816. Find a nominal voltage that would give luminance L. Failing that, it returns the
  817. closest possible value. Solves the rectified power law:
  818.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  819.     Rectify(x)=x if x≥0
  820.     Rectify(x)=0 if x<0
  821. LToVPower takes about 300 microseconds on a Mac II with 8881 arithmetic. This
  822. routine will be called essentially once per clut entry, so computing a whole new
  823. clut with entries 0..255 will take 256*0.3 ms = 77 ms. Some experiments use only
  824. part of the clut, and will therefore take less time. In experiments that do use
  825. the whole clut it may be desirable to compute the clut table in advance.
  826. dgp 10/28/89
  827. */
  828. {
  829.     double V;
  830.     
  831.     if(L<LP->power[0]) L=LP->power[0];
  832.     V=(pow(L-LP->power[0],1.0/LP->power[3])-LP->power[1])/LP->power[2];
  833.     if(V>LP->VMax) V=LP->VMax;
  834.     return V;
  835. }
  836.  
  837. double LToVQuadratic(LuminanceRecord *LP,double L)
  838. /*
  839. This is a quick, approximate version of LToV() that is used by LToVPolynomial
  840. for an initial guess if no power law fit is available. Find nominal voltage that
  841. would give luminance L. Since a quadratic fit to the luminance calibration is
  842. only a fair approximation, this routine serves only to find a quick approximate
  843. answer to the root of a higher order polynomial fit. Solves the quadratic
  844. equation by method recommended in Numerical Recipes in C. My choice of root
  845. assumes that L is an increasing function of V over the operating range of the
  846. display. If concave (i.e. q[2]<0) use the larger root. If convex (i.e. q[2]>0)
  847. use the smaller root. In plain English, if the parabola is right side up then we
  848. take the right hand branch. If the parabola is upside down then we take the left
  849. hand branch. dgp 8/12/89
  850. */
  851. {
  852.     register double a,b,c,q;
  853.     double V;
  854.  
  855.     if(L>LP->LMax) return LP->VMax;
  856.     if(L<LP->LMin) return LP->VMin;
  857.     c=LP->q[0]-L;
  858.     b=LP->q[1];
  859.     a=LP->q[2];
  860.     if(a == 0.0) return -c/b;
  861.     if(b<0.0)
  862.         q=-0.5*(b-sqrt(b*b-4.0*a*c));
  863.     else
  864.         q=-0.5*(b+sqrt(b*b-4.0*a*c));
  865.     if(a*b<0.0)
  866.         V=q/a;
  867.     else
  868.         V=c/q;
  869.     return V;
  870. }
  871.  
  872. double SetLuminanceRange(LuminanceRecord *LPtr,double lowLuminance,double highLuminance)
  873. /*
  874. The user will never call this routine directly, as it is called automatically by
  875. SetLuminance and SetLuminances. This routine does the hard work of figuring out
  876. which DACs to fix to optimally represent a given luminance range. The goal is to
  877. fix the coarser DACs and vary the finer DACs to represent the image. The CHANGES
  878. in luminance (i.e. the contrast) will have the accuracy of the coarsest of the
  879. variable DACs. To avoid the overhead of figuring all this out every single time
  880. you call SetLuminance(), the results are stored temporarily in your
  881. LuminanceRecord, and are only recomputed if you request a different luminance
  882. range.
  883.  
  884. This function checks and returns immediately if you call it with the same
  885. luminance range that it already has stored in your LuminanceRecord.
  886. */
  887. {
  888.     register LuminanceRecord *LP=LPtr;
  889.     register int i;
  890.     double dL,dL2,LOffset;
  891.     double VSum,VFixed;
  892.     int dacs,fixed,dac[DACS];
  893.     double gain[DACS],V[DACS];
  894.     double VGoal;
  895.     double lowV,highV;
  896.     double tolerance;
  897.     double VHalfStep;
  898.     double dV;
  899.     unsigned short *RGBArray;
  900.     
  901.     /* -1. If necessary, swap low & high */
  902.     if(lowLuminance > highLuminance){
  903.         dL=lowLuminance;
  904.         lowLuminance=highLuminance;
  905.         highLuminance=dL;
  906.     }
  907.  
  908.     /* 0. Return right away if our work has already been done */
  909.     if(LP->rangeSet == luminanceSet &&
  910.         LP->lowLuminance==lowLuminance &&
  911.         LP->highLuminance==highLuminance)
  912.             return LP->tolerance;
  913.  
  914.     /* 1. sort the DACs by gain, discard any with zero or negative gain.
  915.     Gain is defined as the change in V when a single DAC is increased from 0 to 1.
  916.     Sort the gains so that gain[0]≥gain[1]≥gain[2]...
  917.     Note this order is opposite to that used in Pelli & Zhang (1991). The answers
  918.     are the same, but the notation is different.
  919.     */
  920.     if(LP->rangeSet != luminanceSet){
  921.         /* This only needs to be done once, even if the range changes. */
  922.         gain[0]=LP->r;
  923.         gain[1]=LP->g;
  924.         gain[2]=LP->b;
  925.         #if DACS>3
  926.             #error "Current implementation doesn't support more than 3 DACs."
  927.         #endif
  928.         for(i=0;i<DACS;i++)dac[i]=i;
  929.         Sort(DACS,gain,dac);            /* sort so that gain[0]≥gain[1]≥gain[2] */
  930.         dacs=0;
  931.         for(i=0;i<DACS;i++) if(gain[i]>0.0)dacs++;
  932.         VHalfStep=gain[dacs-1]/2.0;        /* half the smallest step */
  933.         for(i=0;i<COLORS;i++)LP->table[i].value=i;
  934.         LP->dacSize=Log2L(LP->VMax*2-1);    // in bits, rounded up
  935.         if(LP->dacSize!=8)printf(
  936.             "Caution: Your video card seems to have %ld-bit video dacs, and Luminance.c\n"
  937.             "has only been tested with 8-bit dacs.\n",LP->dacSize);
  938.         /* Apple's convention is to provide a 16-bit value to the dac. Luminance.c
  939.         saves time and space by computing values with only as much precision as
  940.         needed for the the video card's dac (i.e. 0 to VMax, or dacSize bits). This
  941.         value must be shifted left by leftShift bits to yield a standard 16-bit value.
  942.         Note that Apple normally recommends scaling the value to fill the entire
  943.         16-bit range, e.g. multiplying an 8-bit value by 0x101, but that doesn't
  944.         seem appropriate here, where we know the dac size, and are more concerned with
  945.         step size that with range.
  946.         */
  947.         LP->leftShift=16-LP->dacSize;
  948.     }
  949.     else {
  950.         for(i=0;i<DACS;i++){
  951.             gain[i]=LP->gain[i];
  952.             dac[i]=LP->dac[i];
  953.         }
  954.         dacs=LP->dacs;
  955.         VHalfStep=LP->VHalfStep;
  956.     }
  957.     for(i=0;i<DACS;i++)V[i]=0.0;    /* initialize fixed dac voltages */
  958.     
  959.     /* 2. Transform lowLuminance and highLuminance to lowV and highV. */
  960.     lowV=LToV(LP,lowLuminance);        /* takes 0.15 ms */
  961.     highV=LToV(LP,highLuminance);    /* takes 0.15 ms */
  962.  
  963.     /* 3. Designate the finest dacs as variable until we have enough to cover
  964.     the range, then designate the rest as fixed. (DACs 0..fixed-1 will be fixed.)
  965.     */
  966.     VGoal=fabs(highV-lowV);
  967.     VSum=0.0;
  968.     fixed=0;
  969.     for(i=dacs-1;i>=0;i--) {
  970.         if(VSum >= VGoal)fixed++;
  971.         VSum += gain[i]*(LP->VMax - LP->VMin);
  972.     }
  973.     /* sum gains of all variable dacs. This is called gVary by Pelli & Zhang (1991) */
  974.     tolerance=0.0;
  975.     for(i=dacs-1;i>=fixed;i--)tolerance+=gain[i];
  976.     /* scale by highest luminance gain in the requested range */
  977.     dL=fabs(VToL(LP,lowV+1.0)-VToL(LP,lowV));
  978.     dL2=fabs(VToL(LP,highV)-VToL(LP,highV-1.0));
  979.     tolerance *= MAX(dL,dL2);
  980.  
  981.     /* 4. Temporarily set the variable DACs to their midpoints. Now set the fixed DACs
  982.     to most accurately represent (lowV+highV)/2.
  983.     */
  984.     VGoal=(lowV+highV)/2.0;
  985.     VSum=0.0;
  986.     VFixed=(LP->VMax + LP->VMin)/2.0;
  987.     for(i=fixed;i<dacs;i++) VSum+=gain[i]*VFixed;    /* The mid-voltage of the variable DACs */
  988.     for(i=0;i<fixed;i++) {
  989.         V[i]=floor((VHalfStep+VGoal-VSum)/gain[i]);    /* the unattenuated voltage of DAC i */
  990.         V[i]=MAX(LP->VMin,MIN(LP->VMax,V[i]));
  991.         VSum += V[i]*gain[i];
  992.     }
  993.     VSum=0.0;
  994.     for(i=0;i<fixed;i++) VSum+=gain[i]*V[i];
  995.     VFixed=VSum;
  996.     
  997.     /* 5. The limited precision of the fixed DACs may result in a 
  998.     small offset between the requested and now available range. The offset will be
  999.     at most half a step of the finest of the fixed DACs. It seems reasonable
  1000.     to offset the requested luminance range by up to that small amount to better fit
  1001.     the now available range. */
  1002.     LOffset=0.0;
  1003.     if(fixed>0){
  1004.         VSum=VFixed;
  1005.         for(i=fixed;i<dacs;i++) VSum += LP->VMin*gain[i];
  1006.         dV= VSum - MIN(lowV,highV);
  1007.         dV= MIN(dV,gain[fixed-1]/2.0);
  1008.         if(dV>0.0) LOffset+=VToL(LP,VSum+dV)-VToL(LP,VSum);
  1009.         VSum=VFixed;
  1010.         for(i=fixed;i<dacs;i++) VSum += LP->VMax*gain[i];
  1011.         dV= MAX(lowV,highV) - VSum;
  1012.         dV= MIN(dV,gain[fixed-1]/2.0);
  1013.         if(dV>0.0) LOffset-=VToL(LP,VSum+dV)-VToL(LP,VSum);
  1014.     }
  1015.     /* Now save this information in the LuminanceRecord for future use */
  1016.     for(i=0;i<DACS;i++){
  1017.         LP->dac[i]=dac[i];
  1018.         LP->gain[i]=gain[i];
  1019.     }
  1020.     LP->rangeSet=luminanceSet;
  1021.     LP->lowLuminance=lowLuminance;
  1022.     LP->highLuminance=highLuminance;
  1023.     LP->dacs=dacs;
  1024.     LP->VHalfStep=VHalfStep;
  1025.     LP->fixed=fixed;
  1026.     LP->VFixed=VFixed;
  1027.     LP->tolerance=tolerance;
  1028.     LP->LOffset=LOffset;
  1029.  
  1030.     /* Copy into FIXED variables */
  1031.     for(i=0;i<DACS;i++)LP->_gain[i]=DoubleToFix(LP->gain[i]);
  1032.     LP->_VHalfStep=DoubleToFix(VHalfStep);
  1033.     LP->_VFixed=DoubleToFix(VFixed);
  1034.     LP->_tolerance=DoubleToFix(tolerance);
  1035.     LP->_LOffset=DoubleToFix(LOffset);
  1036.  
  1037.     /* cache the fixed DAC values, and zero the rest for good measure */
  1038.     RGBArray = (unsigned short *) &LP->rgb;    /* treat as an array */
  1039.     for(i=0;i<DACS;i++){
  1040.         RGBArray[LP->dac[i]] = (short)V[i]<<LP->leftShift;
  1041.     }
  1042.     return tolerance;
  1043. }
  1044.  
  1045. static void Sort(int n,double arr[],int krr[])
  1046. /*
  1047. Slightly modified sort routine piksr2() from Numerical Recipes in C. I changed
  1048. "float" to "double", and I changed second array to int. I reversed the order, so
  1049. largest element would be first. I changed it to use conventional c arrays,
  1050. starting at index 0.
  1051. */
  1052. {
  1053.     register int i,j;
  1054.     double a;
  1055.     int k;
  1056.  
  1057.     for(j=1;j<n;j++) {
  1058.         a=arr[j];
  1059.         k=krr[j];
  1060.         i=j-1;
  1061.         while (i >= 0 && arr[i] < a) {
  1062.             arr[i+1]=arr[i];
  1063.             krr[i+1]=krr[i];
  1064.             i--;
  1065.         }
  1066.         arr[i+1]=a;
  1067.         krr[i+1]=k;
  1068.     }
  1069. }
  1070.